/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/************************************************************************************************************
 * lpp.c -- Direct evaluation of the intramembrane pressure ("lateral pressure profile") 
 ************************************************************************************************************/

#include "dataStructs.h"

//#define SPLIT_LATERAL_NORMAL_COMPONENTS // uncomment to calculate the individual components (normal and lateral)

extern const VecR *siteCoords; 
extern const VecR *siteVelocity; 
extern const Site *site; 
extern const VecR region; 
extern const int nSites, stepCount, sizeHistLpp, limitLpp, nDOPCsDSPCs, nDOPEs,  
  halfSizeHistLpp; // defined in startUp.c
extern const real deltaT, *mass, timeNow, lppDeltaZ, extTemperature;
extern int countLpp; // global, defined & initialised in startUp.c
extern real *totalHistLpp, *histLppP1, *histLppP2, *histLppNormalComponent, *histLppLateralComponent,
  meanLppSlabVol, meanLppDeltaZ, meanLppBoxHeight; // global, defined and initialised in startUp.c

void ZeroLppHistograms(); // external function declaration (function defined in startUp.c)

/************************************************************************************************************
 * DistributeVirial -- local stress decomposition
 * The Harasima contour is normally used -- "#define IRVING_KIRKWOOD" is normally commented out            
 * To use the Irving-Kirkwood contour, restore "#define IRVING_KIRKWOOD"
 ************************************************************************************************************/
void DistributeVirial( const RMat virTens, const int i, const int j ) 
{       
  int iSlab, jSlab;
  real latPress = ExtractLatPress( virTens.u ); // pL - pN

  //#define IRVING_KIRKWOOD // uncomment to use IK contour (default is Harasima)
#ifdef IRVING_KIRKWOOD
  int rightmostSlab, leftmostSlab, nInterSlabs, slabCounter;
  real zij, fLatPress;
  VecR rij;
  VSub( rij, site[i].r, site[j].r );
  VWrap( rij, z );   /* Minimum Image Convention */
  zij = fabs( rij.z );
#endif

#ifdef SPLIT_LATERAL_NORMAL_COMPONENTS  
  real normPressComp = ExtractNormPressComponent( virTens.u ); // pN
  real latPressComp = ExtractLatPressComponent( virTens.u ); // pL
#endif

  iSlab = siteCoords[ i ].z / lppDeltaZ + halfSizeHistLpp; // make all slab bins positive
  jSlab = siteCoords[ j ].z / lppDeltaZ + halfSizeHistLpp; 

#ifdef IRVING_KIRKWOOD
  if ( iSlab == jSlab ) /* if both sites are in the same slab */
    totalHistLpp[ iSlab ] += latPress; /* place the entire virial in the common slab */
  else {
    rightmostSlab = Max( iSlab, jSlab );
    leftmostSlab = Min( iSlab, jSlab ); 
    nInterSlabs = ( rightmostSlab - leftmostSlab ) - 1; // compute the #slabs between the sites, without considering MIC/PBC
    if ( nInterSlabs == 0 || nInterSlabs == (sizeHistLpp-2) ) { // if directly adjacent, or through PBC
      fLatPress = .5 * latPress; // assign half virial to each of the two slabs
      totalHistLpp[ iSlab ] += fLatPress; 
      totalHistLpp[ jSlab ] += fLatPress; 
    } else { // this is the general case when the two sites are separated by one or more slabs
      if ( nInterSlabs < halfSizeHistLpp ) { // if PBC are definitely not involved  
	fLatPress = latPress / (nInterSlabs+2);  // then, distribute this to the inter-slabs
	totalHistLpp[ iSlab ] += fLatPress; // in primis, assign the adequate bit of virial to iSlab and jSlab
	totalHistLpp[ jSlab ] += fLatPress; 
	for ( slabCounter = leftmostSlab + 1; slabCounter < rightmostSlab; slabCounter++ ) // loop through the inter-slabs between i and j
	  totalHistLpp[ slabCounter ] += fLatPress;
      } else { // PBC wrapping 
	fLatPress = latPress / (sizeHistLpp-nInterSlabs);
	totalHistLpp[ iSlab ] += fLatPress; // in primis, assign the adequate bit of virial to iSlab and jSlab
	totalHistLpp[ jSlab ] += fLatPress; 
	for ( slabCounter = 0; slabCounter < leftmostSlab; slabCounter++ ) // scan the bottom slabs 
	  totalHistLpp[ slabCounter ] += fLatPress;
	for ( slabCounter = rightmostSlab + 1; slabCounter < sizeHistLpp; slabCounter++ ) // scan the top slabs
	  totalHistLpp[ slabCounter ] += fLatPress;
      }
    } 
  } 
#else /* Harasima contour */
  totalHistLpp[ iSlab ] += .5 * latPress;
  totalHistLpp[ jSlab ] += .5 * latPress;
#endif
#ifdef SPLIT_LATERAL_NORMAL_COMPONENTS  
  histLppNormalComponent[ iSlab ] += .5 * normPressComp;
  histLppNormalComponent[ jSlab ] += .5 * normPressComp;
  histLppLateralComponent[ iSlab ] += .5 * latPressComp;
  histLppLateralComponent[ jSlab ] += .5 * latPressComp;
#endif
}

static void PrintLpp()
{
  int bin;
  real distanceFromBilayerCentre;
  FILE *lppFilePtr;
  char lppFileName[ 20 ];
  char time_ns[ 9 ]; // string containing timeNow, e.g. "100", in ns

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "lpp-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( totalHistLpp[ bin ] * PRESSURE_atm ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "P1-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppP1[ bin ] * PRESSURE_atm * LENGTH_nm ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "P2-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.; 
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppP2[ bin ] * PRESSURE_atm * Sqr( LENGTH_nm ) ) ); 
    fprintf( lppFilePtr, "\n" ); 
  }
  fclose( lppFilePtr );

#ifdef SPLIT_LATERAL_NORMAL_COMPONENTS  
  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "PN-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.;
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppNormalComponent[ bin ] * PRESSURE_atm ) );
    fprintf( lppFilePtr, "\n" );
  }
  fclose( lppFilePtr );

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  lppFileName[ 0 ] = '\0';
  strcat( lppFileName, "PL-" );
  strcat( lppFileName, time_ns );
  strcat( lppFileName, ".dat" );
  lppFilePtr = fopen( lppFileName, "w" );
  
  for ( bin = 0; bin < sizeHistLpp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanLppDeltaZ - meanLppBoxHeight / 2.;
    fprintf( lppFilePtr, "%7.3f", ( ( distanceFromBilayerCentre ) * LENGTH_nm ) );
    fprintf( lppFilePtr, "%15.3f", ( histLppLateralComponent[ bin ] * PRESSURE_atm ) );
    fprintf( lppFilePtr, "\n" );
  }
  fclose( lppFilePtr );
#endif
}

static void ComputeSurfaceTension() { // integrate the lateral pressure profile: gamma = - 0.5 * int( lpp ) dz
  int bin;
  FILE *gammaFilePtr;
  real gammaLeft = 0.;  // surface tension of "left" monolayer
  real gammaRight = 0.; // surface tension of "right" monolayer
  
  for ( bin = halfSizeHistLpp - 1; bin >= 0; bin-- )  // "left" monolayer - going backwards from bilayer centre to water
    gammaLeft  -= totalHistLpp[ bin ];   // accumulate pN-pL = -(pL-pN) = -pi
  
  for ( bin = halfSizeHistLpp; bin < sizeHistLpp; bin++ ) // "right" monolayer
    gammaRight -= totalHistLpp[ bin ];   // accumulate pN-pL = -(pL-pN) = -pi
    
  gammaLeft *= meanLppDeltaZ; // normalise
  gammaRight *= meanLppDeltaZ; // normalise

  gammaFilePtr = fopen( "surfaceTensionPerMonolayer.log", "a");
  fprintf( gammaFilePtr, "Surface tension of \"left\"  monolayer from integrating lpp: gammaLeft  = %10.6f dyn/cm\n",   ( gammaLeft  * TENSION_dyn_cm ) );
  fprintf( gammaFilePtr, "Surface tension of \"right\" monolayer from integrating lpp: gammaRight = %10.6f dyn/cm\n\n", ( gammaRight * TENSION_dyn_cm ) );
  fclose( gammaFilePtr );
} 

/********************************************************************************************************************************************
 * ComputeLatPressMoments -- computes the first & second moment of the lateral pressure profile [Cantor, Biophys J, 80, 2284 (2001)]
 ********************************************************************************************************************************************/
static void ComputeLatPressMoments() 
{
  const real T_room_K = 293;  // room temperature [K], same as 20 Celsius
  const real kBT_room_J = kB_IN_J_OVER_K * T_room_K; // thermal energy at room temperature ~ 4.0453E-21 J

  const real DOPC_BILAYER_BENDING_MODULUS_J = 0.76e-19;  // [J] from Pan et al., BJ [2008]
  const real DOPE_BILAYER_BENDING_MODULUS_J = 0.90e-19;  // [J] from Chen & Rand, BJ [1997]
  /*   const real c0_DOPC_per_nm = -0.115; // [1/nm] DOPC spontaneous curvature, from Chen & Rand, BJ, [1997] */
  /*   const real P1_DOPC_kBT_per_nm = 0.5*DOPC_BILAYER_BENDING_MODULUS_J/kBT_room_J*c0_DOPC_per_nm; // pressure profile's first moment */

  real p, P1, P2, kG_left, kG_right, kG, P1Left, P1Right, P2Left, P2Right, z;
  real c0_left, c0_right, c0; // monolayer spontaneous curvatures: left monolayer, right monolayer, average

  real xi = 1.; // [nm] z-distance of the pivotal plane (polar-apolar interface) from the bilayer centre [Orsi et al, JPCM, 2010]
  real kBendingMonolayer_J; /* monolayer bending rigidity modulus in units of J */

  int bin;
  FILE *fPtr;
  
  P1Left = P1Right = P2Left = P2Right = 0.;

  // Setting the *monolayer* bending modulus for subsequent c0 calculation - this formula gives a rough estimate for mixtures too
  kBendingMonolayer_J = 0.5 * ( nDOPCsDSPCs * DOPC_BILAYER_BENDING_MODULUS_J + nDOPEs * DOPE_BILAYER_BENDING_MODULUS_J ) / ( nDOPCsDSPCs + nDOPEs );
  //  printf ("***Just checking - kBendingMonolayer_J = %e\n",kBendingMonolayer_J );
  
  for ( bin = halfSizeHistLpp - 1; bin >= 0; bin-- ) { // "left" monolayer - going backwards from bilayer centre to water
    z = ( ( halfSizeHistLpp - 1 ) - bin + 0.5 ) * meanLppDeltaZ;  // compute the (positive!) distance "z" from the bilayer centre
    //    printf("\nz = %f\n", z * LENGTH_nm);
    p = totalHistLpp[ bin ];   // get the local net lateral pressure acting at distance "z"
    histLppP1[ bin ] = - z * p; 
    P1Left -= histLppP1[ bin ]; //accumulate first integral moment of the lateral pressure (negative summation as we are on the negative side of the z-axis)
    histLppP2[ bin ] = Sqr( z ) * p;
    P2Left += histLppP2[ bin ]; // accumulate second integral moment of the lateral pressure
  }
  
  P1Left *= meanLppDeltaZ;
  P2Left *= meanLppDeltaZ;

  c0_left = P1Left / ( kBendingMonolayer_J / ENERGY_J ) / LENGTH_nm; // units: [ 1 / nm ] 
  kG_left = 2. * xi * P1Left - P2Left;
  
  fPtr = fopen( "lppMoments.log", "a");
  /*   fprintf( fPtr, "P1 left layer = % 7.4f kB * T / nm\t", P1Left * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "P2 left layer = % 7.4f kB * T\n",     P2Left * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "c0_left = P1_left / k = % 5.4f / nm\t", c0_left ); */
  /*   fprintf( fPtr, "kG_left = 2*xi*P1 - P2_left = % 7.4f kBT\n\n", kG_left * ENERGY_J / kBT_room_J ); */

  for ( bin = halfSizeHistLpp; bin < sizeHistLpp; bin++ ) { // "right" monolayer
    z = ( bin - halfSizeHistLpp + 0.5 ) * meanLppDeltaZ;  // compute the distance from the bilayer centre
    p = totalHistLpp[ bin ];   // get the local net lateral pressure acting at distance "z"
    histLppP1[ bin ] = z * p;
    P1Right += histLppP1[ bin ];
    histLppP2[ bin ] = Sqr( z ) * p;
    P2Right += histLppP2[ bin ];
  }

  P1Right *= meanLppDeltaZ;
  P2Right *= meanLppDeltaZ;
 
  c0_right = P1Right / ( kBendingMonolayer_J / ENERGY_J ) / LENGTH_nm; // units: [ 1 / nm ] 
  kG_right = 2. * xi * P1Right - P2Right;

  /*   fprintf( fPtr, "P1 right layer = % 7.4f kB * T / nm\t", P1Right * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "P2 right layer = % 7.4f kB * T\n",     P2Right * ENERGY_J / kBT_room_J ); */
  /*   fprintf( fPtr, "c0_right = P1_right / k = % 5.4f / nm\t", c0_right ); */
  /*   fprintf( fPtr, "kG_right = - ( P2_right - 2*xi*P1_right ) = % 7.4f kBT\n\n", kG_right * ENERGY_J / kBT_room_J ); */

  P1 = 0.5 * ( P1Right + P1Left );
  P2 = 0.5 * ( P2Right + P2Left );
  c0 = P1 / ( kBendingMonolayer_J / ENERGY_J ); // units: [ 1 / nm ] 
  kG = 2. * xi * P1 - P2;

  fprintf( fPtr, "Monolayer torque tension tau = P1 = % .2f kB * T_room / nm", P1 * ENERGY_J / kBT_room_J );
  if ( nDOPCsDSPCs ) fprintf( fPtr, " - for DOPC this should range from -0.09 to -2.39 kBT/nm" );
  fprintf( fPtr, "\nMonolayer spontaneous curvature c0 = P1 / k = % .2f / nm (assuming k = %.2e J)\n", c0, kBendingMonolayer_J );
  fprintf( fPtr, "Monolayer Gaussian curvature kG=2*xi*P1-P2=%.2f kBT; should be -%.2f kBT < kG < 0 ", kG*ENERGY_J/kBT_room_J, kBendingMonolayer_J/kBT_room_J ); 
  fprintf( fPtr, "(assuming xi = %.2f nm)\n", xi ); 
  fclose( fPtr );

  fPtr = fopen( "P1.dat", "a");
  fprintf( fPtr, "% 8.6f\n", P1 * ENERGY_J / kBT_room_J );  // kB * T_room / nm
  fclose( fPtr );

  fPtr = fopen( "P2.dat", "a");
  fprintf( fPtr, "% 8.6f\n", P2 * ENERGY_J / kBT_room_J ); //  kB * T
  fclose( fPtr );

} 

void EvalLpp() 
{
  RMat kinEnDensTensor; // kinetic energy tensor density tensor - used to in the evaluation (accumulation) of the pressure profile kinetic contribution 
  int n, bin; // "bin" is the histogram bin number

  real slabVolume = region.x * region.y * lppDeltaZ;   // computing the volume over which the density is worked out
  
  DO_SITE { // computing and sorting kin en contribution to (pL - pN)
    VOuter( kinEnDensTensor.u, siteVelocity[ n ], siteVelocity[ n ] );
    MScale( kinEnDensTensor.u, mass[ site[ n ].type ] );
    bin = siteCoords[ n ].z / lppDeltaZ + halfSizeHistLpp;    
    totalHistLpp[ bin ] += ExtractLatPress( kinEnDensTensor.u );
#ifdef SPLIT_LATERAL_NORMAL_COMPONENTS  
    histLppNormalComponent[ bin ]+= ExtractNormPressComponent( kinEnDensTensor.u );
    histLppLateralComponent[ bin ]+= ExtractLatPressComponent( kinEnDensTensor.u );
#endif
  }
  
  meanLppSlabVol += slabVolume; // accumulating - really only necessary for NPT (redundant at constant volume, but its no big deal 2 compute, hence leave it)
  meanLppDeltaZ += lppDeltaZ;  // lppDeltaZ computed in SingleStep()
  meanLppBoxHeight += region.z; 
  
  ++countLpp;

  if ( limitLpp == countLpp ) {  // if happy with the amount of data collected
    for ( bin = 0; bin < sizeHistLpp; bin++ )  {// normalising to have an averaged value
      totalHistLpp[ bin ] /= meanLppSlabVol; 
#ifdef SPLIT_LATERAL_NORMAL_COMPONENTS  
      histLppNormalComponent[ bin ]  /= meanLppSlabVol; 
      histLppLateralComponent[ bin ] /= meanLppSlabVol; 
#endif
    }
    
    meanLppDeltaZ /= countLpp;
    meanLppBoxHeight /= countLpp; 
    
    ComputeSurfaceTension();
    ComputeLatPressMoments();

    PrintLpp();

    countLpp = 0;
    meanLppSlabVol = 0.;
    ZeroLppHistograms();
  }
}
